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Abstract. We present a semi-analytic approach to solv- 
ing the Boltzmann equation describing the comptonisa- 
tion of low frequency input photons by a thermal distri- 
bution of electrons in the Thomson limit. Our work is 
based on the formulation of the problem by Titarchuk 
& Lyubarskij ( 1995Q , but extends their treatment by ac- 
commodating an arbitrary anisotropy of the source func- 
tion. To achieve this, we expand the cigcnfunctions of the 
integro/differential eigenvalue problem defining the spec- 
tral index of comptonised radiation in terms of Legendre 
polynomials and Chebyshev polynomials. The resulting 
algebraic eigenvalue problem is then solved by numerical 
means, yielding the spectral index and the full angular 
and spatial dependence of the specific intensity of radi- 
ation. For a thin (to < 1) plasma disk, the radiation is 
strongly collimated along the disk surface - for an optical 
thickness of To = 0.05, the radiation intensity along the 
surface is roughly ten times that along the direction of the 
normal, and varies only slightly with the electron temper- 
ature. Our results for the spectral index confirm those of 
Titarchuk & Lyubarskij over a wide range of electron tem- 
perature and optical depth; the largest difference we find 
is roughly 10% and occurs at low optical depth. 

Key words: Radiative transfer - Radiation mechanisms: 
miscellaneous - Methods: analytical - Galaxies: active - 
X-rays: galaxies - X-rays: stars 



1. Introduction 

The process of Comptonisation, in which soft photons in- 
crease their energy by scattering in a gas of hot electrons is 
thought to be the main mechanism responsible for the for- 
mation of non-thermal continuum spectrum in a range of 
astrophysical objects. Much of the early work on this sub- 
ject assumed that the optical depth of the scattering cloud 



was fairly large, and that the photon frequency v and elec- 
tron temperature T c were both small: x = hv/m c c 2 <C 1, 
O = k-QT e /m c c 2 <C 1. The transport of a photon in both 
configuration space and in energy space can then be ap- 
proximated by a Fokker-Planck equation and the compu- 
tation of the spectrum can be split into two parts - the 
evaluation of the distribution of the number of scatterings 
undergone by a photon, and the convolution of this with 
the solution for the evolution of the spectrum in an infinite 



homogeneous medium (Sunyaev & Titarchuk 1980). 

In the wake of this important result, many authors 
have investigated methods of extending the range of ap- 
plicability of the calculations, e.g., by calculating correc- 
tions to the diffusion coefficient in energy space to higher 
order in the parameters 9 and x (Prasad et al. 1988; see 
also Cooper 1971). Of particular interest is the generali- 
sation to scattering media of optical depth r ~ 1, since 
such conditions are indicated in many applications (e.g., 
AGN: Haardt et al. [1994| Zdziarski et al. |l995|an d galactic 
black hole candidates Sunyaev & Triimper 1979^ Ebisawa 
et al. [19961 ). Thus, Sunyaev & Titarchuk ( |1985| ) presented 
a numerical solution to the spatial transport of photons 
in slab geometry. Assuming the problem remains separa- 
ble, this can be combined with the solution for diffusion 
in energy space to yield the emergent spectrum. These 
and other generalisations are summarised by Titarchuk 
( 1994 ), who also pointed out that the angular distribu- 
tion of comptonised radiation emerging from an optically 
thin disk forms a 'knife-blade' pattern, collimated almost 
parallel to the surface of the disk. In addition to analytic 
work, comptonisation has been investigated using numer- 
ical methods (Katz 1976 ; Poutanen & Svensson |l996| ), in 
particular the Monte-Carlo simulation technique (Pozd- 
nyakov et al. |1983|; Zdziarski [198$ Hua & Titarchuk [1995; 
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Stern et al. [L995a| , [I995b|) . 

In a recent paper, Titarchuk & Lyubarskij ( 1995 ) have 
attacked the problem using the Boltzmann equation, with- 
out recourse to a Fokker-Planck approximation in either 
configuration or energy space. As well as confirming the 
result that power-law spectra are produced when low fre- 
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qucncy photons are scattered in the Thomson regime (i.e, 
when the dimensionless photon energy x' measured in the 
electron rest frame satisfies x' <C 1), they demonstrate ex- 
plicitly the separation of the problem into its configuration 
and energy space parts, provided that the source junction 
can be considered isotropic. In this case, the problem re- 
duces to an eigenvalue equation for the source function. 
Titarchuk & Lyubarskij also provide analytic approxima- 
tions for the computation of the power-law index a over 
a very large range of optical depth and temperature. 
In this paper we generalise the approach of Titarchuk 



& Lyubarskij (1995) by solving the Boltzmann equation 
without assuming isotropy of either the radiation or the 
source function. The approach we adopt is to formulate 
the equation determining the power-law index as an inte- 
gral eigenvalue problem. To find the eigenfunctions, we ex- 
pand the angular dependence in a series of Legendre poly- 
nomials, and the spatial dependence in a series of Cheby- 
shev polynomials. Our results confirm the accu racy o f the 
formulae presented by Titarchuk & Lyubarskij (1995) and 
also give the full spatial and angular dependences of th e 
comptonised radiation. As predicted by Titarchuk ( 1994 ), 
the radiation from an optically thin slab is strongly colli- 
mated along the surface of the slab. 

The paper is laid out as follows: in Sect. || we present 
the equations leading to the formulation of the integral 
eigenvalue problem for the power-law index a. Section [3] 
(and in more detail the appendix) presents the method 
of solution. The phase function is first of all expanded in 
terms of Legendre polynomials of the scattering angle and 
the integration over electron velocity is performed in the 
case of low temperature (8 C 1) and high temperature 
(O 3> 1). Then, using an expansion in Legendre polyno- 
mials which are complete over the half-range < fi < 1 
of the cosine of the angle between the photon direction 
and the normal to the slab, the problem is converted to a 
system of differential equations in the spatial coordinate 
normal to the disk surface. This is finally reduced to an 
algebraic eigenvalue problem by expanding the spatial de- 
pendence in a series of Chebyshev polynomials. Singular 
value decomposition is then used to find the eigenvalues 
and eigenfunctions. Our results are presented in Sect. |j. 
These consist of plots of the spectral index as a function of 
temperature O and half-thickness tq of the slab, the angu- 
lar dependence of the specific intensity of radiation, and 
its spatial distribution. We compare these to the formulae 
for a presented by Titarchuk & Lyubarskij (1995). Sec- 
tion |5| contains a summary of our conclusions, and a short 
discussion of their range of applicability in astrophysical 
sources. 



2. Formulation of the basic equations 

Let us first define the geometry of the problem. We con- 
sider an infinite disk of thickness 2zq containing a uni- 
form, non-degenerate gas of free-electrons of temperature 



T and number density n c . The only process of importance 
for the transport of photons in this disk is Compton scat- 
tering. The optical half-thickness of the disk is defined as 
To = aTn c zo, where or = 6.65 x 10 -25 cm 2 is the Thomson 
cross-section. Let the spatial coordinate normal to the disk 
be z, with z — on the mid-plane of the disk, and define 
the optical depth variable r = n e arz which is bounded by 

-TO < T < Tq. 

We denote the cosine of the angle between the normal 
to the disk and the direction of a photon after scattering 
by ll. In the same sense, the direction before a scattering- 
event is denoted by Let rj be the cosine of the angle 
between these directions. For an isotropic electron distri- 
bution, the phase function, which describes the change in 
photon direction due to scattering, depends only on r\. To 
calculate this phase function, we have to perform an in- 
tegral over the electron velocity. Because of the isotropy 
of the electron distribution, the electron direction does 
not refer to the disk normal. It is therefore convenient 
to switch to a coordinate system, in which the electron 
direction defines the z-axis. The cosines of the photon di- 
rections with respect to the electron direction are then de- 
noted by /x (after scattering), and fii (before scattering). 
An analogous notation is used for the azimuthal angles. 

We are interested in a situation in which low frequency 
radiation is injected into the disk, is scattered by the 
electrons, and forms a power-law spectr um at high fre- 
quency, as shown by Shapiro et al. (1976) and Sunyaev & 
Titarc huk ( 1980 ), and seen in numerical studies of Katz 
( 1976 ). In the power-law regime, there is no source of ra- 
diation in the disk, and no radiation which enters from 
the outside. The time independent, polarisation averaged, 
equation of transfer for the specific (up-scattered) inten- 
sity I(v, fj,, t) is then given by 



9 Tf \ 1 

fJ,—I{U,fJ,,T) = 

OT Il c (JT 



dv± J dili 

n 4tt 



-a s (vi — > v, rf)I(ui, Hx,t) - a s (v — > v\, rj)I(v, ll, t) 



(1) 



see e.g., Pomraning 1973). This transport equation is a 



special case of the linearised Boltzmann equation and we 
shall simply refer to it as the Boltzmann equation. 

The power-law part of the spectrum we describe occurs 
at frequencies lower than that of the Wien cut-off (hv < 
kBT c ), so that the energy change of the photon due to 
the recoil of the electron, can be neglected in comparison 
with the Dopplcr shift of the photon. In the electron rest 
frame (in which quantities are adorned with a prime) , the 
classical Thomson scattering kernel is then given by 

a' s (u' - i/^rf) = ^-n' c a T {l + (V) 2 } 5{y> ~ i/j) . (2) 

The electrons are distributed isotropically. They are de- 
scribed by a relativistic Maxwell distribution, which is de- 
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fined as 



3. Method of solution 



f(v) 



7 5 exp( 



-7/6) 



47reK 2 (l/e) 



(3) 



where K 2 is the modified Bessel function of the second 



To reduce the integral eigenvalue problem described above 
to an algebraic eigenvalue equation, we first expand the 
phase function in a series of Legendre polynomials: 



kind of order two, and 7 = (1 
following we set c 
according to 



'1/2 



(here and in the 



1). This distribution is normalised R(rj) = uQj(a, Q)Pj{ri) 



(12) 



&%f{v) =4tt / dvv^f(v) = 1. 



(4) 



The scattering kernel in the disk system can be calcu- 
lated by performing a Lorentz boost of a' s , multiplying it 
by f(v) and, integrating over v. Then the scattering kernel 
is given by 



a s (u -> v-L,r], 9) 
■{! + (!- 



3 77, e (7T 
167T VX 
1-7] N 2 



d 6 v 



/(«) 



with the definitions 

D := 1 — jiv , 
1 - /2i u . 



r)} 



•5 



7 

. V xi a; 



(5) 



(6) 



We now look for a solution of the Boltzmann equation of 
the form 



J{v,fj,,r) =J{h,t)x a 



(7) 



Inserting this, and Eq. (5), into Eq. (1), and performing 
the trivial integral over v\ using the S- function, we find 



where the source function B{]i, t) is defined as 

1 2rr 

B (^ r ) = ^ J d Mi J d<j}R(rj) J(^i,t) . 

-1 

The phase function R{rf) is given by 



(8) 



(9) 



(10) 



The last three equations (8, 9 and 10), together with the 
boundary condition 



J(/x < 0, r = r ) = = J(fi > 0, r = -to) : 



define an integral eigenvalue problem. Our aim is to reduce 
this to an algebraic eigenvalue equation by expanding the 
intensity J(/z,t) into a polynomial series, as mentioned 
above and shown in the Sects. |A.3| and A. 4. We further 



perform the phase function integral using an expansion 
for -c 1 and 8>1. This is shown in the Sects. A.l and 



1 See also Titarchuk & Lyubarskij (1995) Eq. (2) after cor- 
rection of a minor typographical error. 



i=0 



The details of this and the following steps are shown in 
the appendix. The coefficients u>i(a, O) can be written, 
in the case of small electron temperature (9 <C 1), as 
polynomials in a, with coefficients in terms of moments 
of the Maxwell distribution, which can easily be calcu- 
lated numerically. In the case of high electron temperature 
(6 > 1), the coefficients u>i(a, 9) can be expressed using 
the incomplete T-function. 

The essential step is to expand the specific intensity 
J(/i, t) into a series of Legendre and Chebyshev polyno- 
mials. For the half-space < fj, < 1, we write (compare 
Eq. (A23)) 



J(/x,t)| 



EZn 

n=0 



1 



P„(2 M -l)Q+(r) 



(13) 



We now expand each of these N + 1 expansion coeffi- 
cients Qni T ) m to a series of Chebyshev polynomials (Tj). 
In vector- form this reads (compare Eq. (A45)) 



K 



Q + (t/t ) = q + ^[l-7 1 i (-T/T )] q t 



(14) 



All expansion coefficients can be represented by a common 
(K + 1) • (N + 1) dimensional vector: 



fq \ 

w 



(15) 



As shown in detail in the appendix, the eigenvalue prob- 
lem (Eqs. 8-10) becomes with these expansions a set of 
homogeneous linear equations for the vector q, which is 
very easy to treat numerically. The set of equations can 
be written as a matrix equation (see Eq. (A58)): 



(11) F(a,9,T ) 9 =0. 



(16) 



The first step is to calculate the matrix F(a, 9,Ta) for 
given values of t and 9. The solution for the spectral 
index a can be found by solving the following equation 
for the determinant of F(a) : 



det [£(«)] = 0. 



(17) 



We used Mathematica to find the roots of this equation, 
which give the spectral index a. For an expansion of the 
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to 



0=5/511 




Note that the spectral index is the exponent of the in- 
dependent functions D and D\ in the phase function. For 
non-relativistic plasma temperatures, the spectral index 
can be well above 1. In this case, accuracy is achieved only 
if the phase function is expanded to high order in a Taylor 
series in v — we use an expansion up to 16th. order [L = 16 
in Eq. (A5)). This sum is represented by an expansion to 
10th. order in r] (M = 10 in Eq. (A5)). In the relativistic 
temperature regime, on the other hand, the spectral index 
a is much lower. Therefore, we take into account only the 
leading order of u)i(ot, 0) in the small parameter 1/7. The 
first four coefficients of Eq. (A4) are then given by the 
Eqs. (A17), together with Eq. (A20). The spectral index 



Fig. 1. Spectral index a vs. Thomson optical half thickness 
to, for non-relativistic and relativistic electron plasma temper- 
atures 9 = k B T c /m c « fc B r o /(511keV). 



spatial and angular dependence of J(/x, r) to (e.g.) 9th or- 
der, q is a 100 dimensional vector and F(a) is a 100 x 100 
matrix. A solution of Eq. (17) for this dimension takes less 
than one minute with a Pentium PC in both the relativis- 
tic, and the non-relativistic cases. 

For each value of a, which satisfies Eq. (17), the ex- 
pansion coefficients, and therefore J(/x,r), can be found 
by solving the equation 



£g = 0, 



(18) 



where F is now a known singular square matrix. This ma- 
trix can be decomposed using a si ngular value decompo- 
sition routine (see e.g. Press ct al. 1986| Wolfram 1991 ), 
which yields the vector null-space. This vector can im- 
mediately inserted into Eq. (14), which gives Q ± {t). In- 
serting this into Eq. (13), we end up with the intensity 
J(/x,r). The resulting values for the spectral index a and 
the shape of the intensity in /1 and r are described in 
the next section. Use of the singular value decomposition 
technique has the advantage that an automatic check for 
the presence of repeated roots of Eq. (17) is provided. 
More importantly, however, it immediately provides the 
null-space and the associated specific intensity of radia- 
tion, which must be physically realistic in the sense that 
J(/x, t) must be positive definite. This important condition 
enables one to identify the physically relevant power-law 
index a, which is the smallest positive root of the nonlin- 
ear Eq. (17). 

4. Results 

4-1- Spectral index 

The spectral index a of the comptonised radiation, for 
a given temperature (8 < 1 or 9 > 1), and arbitrary 
optical depth 2tq, can be found by solving Eq. (17). 



„ k B T e 
te> = 

m c 


TO 


M 


a 


50/511 


0.05 





2.83 


50/511 


0.05 


1 


2.80 


50/511 


0.05 


2 


2.75 


50/511 


0.05 


3 


2.74 


50/511 


0.05 


4 


2.74 


50/511 


1.0 





0.678 


50/511 


1.0 


1 


0.656 


50/511 


1.0 


2 


0.652 


50/511 


1.0 


3 


0.652 


50/511 


3.0 





0.186 


50/511 


3.0 


1 


0.179 


50/511 


3.0 


2 


0.179 


4.0 


0.05 





0.366 


4.0 


0.05 


1 


0.361 


4.0 


0.05 


2 


0.359 


4.0 


0.05 


3 


0.359 


4.0 


1.0 





0.0538 


4.0 


1.0 


1 


0.0480 


4.0 


1.0 


2 


0.0480 


4.0 


3.0 





0.0125 


4.0 


3.0 


1 


0.0103 


4.0 


3.0 


2 


0.0103 



Table 1. Spectral index a for two plasma temperatures 
and different expansion parameters M of the source function 
B((j,, r), and for Thomson optical half thickness to = 0.05, 1.0 
and 3.0 . 



is not sensitive to the expansion of the \x and r dependence 
of J(/x,t). It is more than sufficient to choose N = 10 in 
Eq. (A23), and K = 10 in Eq. (A45). Results are given for 
an electron plasma temperature k B T c — 0m c of 5 keV and 
50 keV in the non-relativistic regime. In the relativistic 
regime we choose = 4.0 and = 10.0, corresponding 
roughly to 2 and 5 MeV respectively. Figure [l] shows the 
spectral index a versus the Thomson optical half thickness 
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ro, using the expansion parameters given above (L, M, N 
and K). 

These values of a are in good agreement with those 



given by Titarchuk & Lyubarskij (1995). These authors 
assumed an isotropic source function _B(/i,r) (Eq. 9), 
which is certainly a good approximation for to ^ 1. To re- 
lax this restriction, at least the first three expansion coef- 
ficients of the source function must be taken into account 
(M = 0,1,2), because of the intrinsic 1 + (f/) 2 depen- 
dence of the Thomson scattering kernel. Especially when 
the intensity J(/i, r) is highly anisotropic (which is the 
case for tq C 1, as shown in the next section) the source 
function B(fi,r) depends on the higher expansion coeffi- 
cients of the phase function, which leads to an anisotropy 
of B(h,t). In fact for To = 0.05 the anisotropy of the 
source function at the disk surface for Q = 50/511 be- 
comes 



B(p = 0,t = 0.05) 
B({i = 1,t = 0.05) 



= 2.3 



(19) 



The sensitivity of a to the angular expansion of the source 
function is, however, not very strong, and it is well ap- 
proximated by taking into account the first four Legendre 
polynomials of the expansion (M = 3 in Eq. (A5)). This is 
shown for To = 0.05, 1.0 and 3.0, and two different plasma 
temperatures in Table [l| Note that even for a value of 
a = 2.74 it is necessary to choose a Taylor expansion of 
sufficiently high order (here: L = 16), for any value of M . 



Titarchuk & Lyubarskij ( |1995[ ) suggested interpola- 
tion formulas (see Eq. (17) and Eq. (21) therein, and also 
Eq. (27) in Titarchuk 1994 ), which are valid for all opti- 
cal depths and all plasma temperatures. The values of the 
spectral index given by these formulas, differ by less then 
1% for To » 1 from those given in Fig. |l|; the largest dis- 
crepancy is 10% at smaller values of r (see also Table |l|). 

4-2. Angular distribution 

As described at the end of Sect. ||, a singular value de- 
composition of the matrix F gives J(/x, r) in the form of a 
polynomial of order K in t/tq, and of order N in /z (for 
any set of parameters 6, t and a). We choose L = 16 
and M = 6 (see Eq. (A5)) for the phase function repre- 
sentation, and 9 = 5/511 for the plots discussed below. 
The minimum expansion for the angular and spatial de- 
pendence (K and N) have to be chosen differently for each 
optical depth (see figure captions). 

Because we are interested especially in non-isotropic 
intensities, we define a measure of the anisotropy of 
J(/x, t). First, let J\\ (r) and Jl(t) be the integrated inten- 
sity parallel and perpendicular to the disk, over an interval 
of A// = 0.1: 



o.i 



AM 



d/J,J(fx,r) 



=1 







2t d = 


o.i ; 




T / 


'To = 


1 .o : 




T/ 


't,= 


0.0 1 




Ty 


f r a = 


-1 .o : 













0.1 0.2 0.3 0.4 



1.5 0.6 0.7 0.8 0.9 



Fig. 2. Angular dependence of the specific intensity J(/i, r) for 
a Thomson optical thickness of 2to =0.1. The intensity at the 
surface of the disk is given by t/tq = 1.0. The anisotropy of the 
integrated intensity is A(t ) = {J\\-J±)/{J\\+J±) = 0.82±0.03 
for a wide temperature range. (Expansion up to N — 16 and 
K = 10). 



^ 1 




2t = 


0.4 ; 






't„ = 


1 .o : 


^ 0.6 


/ \ \ T / 


't = 


o.o ; 








-1.0 : 


0.4 








0.2 




















1.5 0.6 0.7 



/J, 



Fig. 3. Angular dependence of the specific intensity J(/i, t) for 
a Thomson optical thickness of 2to = 0.4. The intensity at the 
surface of the disk is given by t/tq = 1.0. The anisotropy of the 
integrated intensity is A(t ) = (J\\-J±)/(J\\+Jx) = 0.47±0.03 
for a wide temperature range. (Expansion up to N — 12 and 
AT = 14). 



J±(t) := J dnJfar) . (20) 

0.9 

From this we define an asymmetryi A, which is for 
isotropic intensity, and 1 for extremely focussed intensity 
along the disk surface: 

J \\( T o) - J±{tq) 



J\\(t ) + Jl(to) 



(21) 



In analogy to the forward-backward asymmetry of the 
electro-weak interaction. 
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Fig. 4. Angular dependence of the specific intensity J(fj., t) for 
a Thomson optical thickness of 2 to = 1.0. The intensity at the 
surface of the disk is given by t/tq = 1.0. The anisotropy of the 
integrated intensity is A(t ) = (J\\-J±)/(J\\+J±) = 0.08±0.02 
for a wide temperature range. (Expansion up to N = 8 and 
K = 16). 



2t = 1.0 




Fig. 5. Spatial distribution of the intensity parallel to the disk 
(as defined in Eq. (20)), normalised to the parallel intensity at 
the disk surface. Note that this plot shows the contributions of 
the half-space < (i < 1. 



At an optical thickness of order unity (2to ~ 1), this 
anisotropy is about zero (see Fig. 0). For smaller values 
of To, A(tq) can become very close to 1, as shown in 
Fig. ^, where J(/i, t) is plotted for an optical depth of 
2tq = 0.1, normalised to the intensity in the middle of 
the disk, parallel to the surface. The boundary condition 
gives J(/i, t = —To) = for < fi < 1 (no radiation en- 
ters the disk from outside, see dashed line). The dotted 
line shows the intensity in the middle of the disk (r = 0), 
whereas the solid line shows the intensity at the surface, 
given by J(/i, r = To). The normalised specific inten- 
sity of the radiation is approximately independent of the 
plasma temperature. As shown by Titarchuk & Lyubarskij 
( |l995|) the electron energy and photon spatial variables 



are completely decoupled if the source function is exactly 
isotropic. However, as discussed in the preceding section, 
the source function has a weak angular dependence, which 
leads to a coupling of the electron energy and photon spa- 
tial variables. This, in turn, implies that the angular de- 
pendence is a function of the plasma temperature. The 
range of variation over the temperature range considered 
is expressed in form of an error of the anisotropy. The 
upper bound of the anisotropy is valid for = 5/511, 
whereas the lower bound was calculated for 9 = 100. 

Even for a moderately small optical depth of 2to = 0.1 
(see Fig.||) the anisotropy is ^4(t ) = 0.82 ± 0.03, which 
means that the radiation in the interval < /i < 0.1 
parallel to the disk surface, is a factor of about 10 more 
intense than the radiation in the interval 0.9 < ji < 1 
perpendicular to the disk. At increased angular resolution 
(smaller A/x) this factor is even larger. 

The reason for the high anisotropy at small optical 
depth is the following: photons which contribute to the 
power-law part of the spectrum have to undergo a num- 
ber of scatterings on electrons to gain the required energy. 
The energy gain has a maximum for back-scattering of 
the photons (A8 = 180°). Thus, those photons most ef- 
fectively boosted in energy and least likely to escape the 
disk are those which move almost parallel to the surface. 
This leads to a strong collimation in the disk plane for an 
optically thin disk. 



4-3. Spatial distribution 

The intensity J(n, r) provides of course not only the an- 
gular distribution, but also the spatial distribution. The 
solution discussed above gives J(fJ>,r) for the half space 
< n < 1. A full angle integrated intensity at any space 
point is then the sum of the two half space intensities, 
which are symmetric with respect to the middle of the 
disk. To compare the spatial distribution to the previous 
figures, we show J(/z, r) also for the half space <u < 1, 
and used the same set of parameters as in Sect. |4.2| . 

Figure ^ shows the intensity J|| parallel to the disk, as 
defined in Eq. (20) for various optical depths, normalised 
to the parallel intensity at the disk surface. At r = —To 
the intensity has to be 0, because we take into account 
only the half space < \x < 1, for which the boundary 
condition is J(/i, r = — t ) = 0. Figure || shows the in- 
tensity Ji perpendicular to the disk, normalised to the 
parallel intensity J|| at the disk surface. Note the differ- 
ence in scale between this plot and Fig. ||. A comparison 
of Fig. H and Fig. |^ provides the anisotropy A(t) for all 
values of r. Because the perpendicular intensity Jj_ drops 
faster for smaller values of t/tq than the parallel intensity 
J||, the anisotropy A(t < tq) is even larger than A(tq), 
which is given in the captions of the Figs. || - [I| 
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Fig. 6. Spatial distribution of the intensity perpendicular to 
the disk (as defined in Eq. (20)), normalised to the parallel 
intensity at the disk surface. Note that this plot shows the 
contributions of the half-space < fj. < 1. 



5. Discussion 

In this paper we have presented a new, semi-analytic 
method of obtaining solutions to the comptonisation prob- 
lem, and used it to find the power-law index of photons 
scattered in the Thomson regime, neglecting the recoil 
of the scattered electron. Because these photons undergo 
a large number of scatterings before emerging from the 
plasma, they have 'forgotten' the details (spatial and an- 
gular dependences) of their injection at low frequency. 
Once in the power-law region, which exists even for op- 
tically thin plasmas, the spatial and angular dependence 
of the specific intensity of radiation ceases to be a func- 
tion of frequency - it is determined by a particular eigen- 
function of the reduced transfer equation and depends 
only on the optical depth and temperature of the plasma. 
We compute this eigenfunction. This distinguishes our 
approach from other numerical techniques in the litera- 
ture. Haardt ( |1993| ), for example, uses an approximation 
method in which anisotropy is accounted for only in the 
first scattering undergone by a photon, so that his results 
are accurate only fairly close to the frequency of injec- 
tion. Poutanen & Svensson (1996) have developed a com- 



prehensive code based on the iterative scattering method, 
which treats the anisotropy 'exactly' (on a discrete grid) 
and accounts for several processes which we neglect. This 
method converges rapidly provided only a few photon 
scatterings are important, i.e., at small optical depth and 
high temperature. However, it would need a large number 
of grid points in order to resolve the sharp dependence 
of the radiation intensity on angle such as displayed in 
Fig. ^. In principle, it should be possible to extend our 
technique to solve 'inhomogeneous' problems, where the 
emergent spectrum depends on the input radiation. How- 
ever, in this case it would be necessary to compute sev- 
eral eigenfunctions. Further extensions of the method to 



geometries other than slab, or to arbitrary electron distri- 
butions are straightforward. 

To apply our results to observations of astrophysical 
objects we have to consider the range of spectral indices 
of order 3 or less, because higher values are probably too 
steep to be observable. A spectral index in this range 
can be achieved by a non-relativistic plasma with optical 
thickness of tq J; 1 , or by a relativistic plasma with optical 
thickness of To ^ 1. The assumption of Thomson scatter- 
ing (as opposed to Klein-Nishina) restricts the relevant 
frequency range to x < 1/(7) (where (7) is the averaged 
Lorentz factor according to Eq. (A6)), i.e. to X-rays for 
the highest temperatures considered here. We also require 
the dominant source of soft photon input to be at a fre- 
quency which is low enough to require more than a single 
scattering before X-ray frequencies are achieved, otherwise 
the pure power-law spectrum is not achieved. The average 
frequency change on scattering is given by {Ax) / x ~ 40 
in the non-relativistic regime, and (Ax)/x ~ (48) 2 for 
relativistic temperatures (see Pozdnyakov et al. 1983 ). 

As an example, consider a plasma with temperature 
fceTc = 50keV and input photons with energy 5eV. After 
20 scatterings, which removes all information of their ini- 
tial distribution, they achieve an energy of roughly 4 keV. 
The resulting normalised photon energy is then x = 0.008, 
which is well below 1/(7) = 0.86. At an optical depth 
of 2to = 0.4 the intensity at the disk surface, parallel 
to it (J||) becomes three times the intensity perpendicu- 
lar to the surface (J_l), {A(r ) = 0.49, compare Fig. |j). 
The spectral index for these values of temperature and 
optical depth is a = 1.77. For relativistic temperatures, 
where spectral indices between and 1 can be produced 
in a thin plasma disk (to <C 1), which leads to a very 
high anisotropy A(to), the input photon energy has to be 
very much lower than 5 eV, in order to have more than 1 
scattering which shifts the photon energy into the X-ray 
regime. 

In particular, our results concerning the degree of 
anisotropy (which is only weakly dependent on temper- 
ature) are relevant to the case of Seyfert galaxies, where 
plasma temperatures of the order of 100 ke V have been 
su ggeste d (Titarchuk & Mastichiadis 1994 ; Zdziarski et 



al. 1995| ). In these objects, the ratio of optical to X-ray 



luminosity should be much smaller for objects seen 'edge- 
on' than for those seen 'face-on' (an effect predicted by 
Haardt & Maraschi 1993 ). It may, however, be difficult to 
disentangle this effect from that of the increased absorp- 
tion of optical radiation expected in edge-on sources. 

Quite apart from application to astrophysically im- 
portant objects, the results of our computations should 
be helpful as a check on other numerical methods of so- 
lution, in particular Monte-Carlo simulations. Here it is 
worth mentioning that we have used a polarisation aver- 
aged treatment of the transport. All current Monte-Carlo 
codes also use this same approximation, so that the re- 
sults they obtain are directly comparable to ours. Gener- 
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alisation to polarisation dependent transfer is in principle 
possible. 
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A. Reduction of the Boltzmann equation to an al- 
gebraic eigenvalue equation 

To express the Boltzmann equation in the form of an al- 
gebraic eigenvalue equation, we express the angular and 
spatial dependence of the source function R(rf) and the 
specific intensity J(/z, r) in terms of Legendre and Cheby- 
shev polynomials. 

Let us first write the source function as a series 



M 



R(rf) = > w,(a, Q)PM, 



i=0 



(Al) 



where the angular dependence of R{if) is expressed in form 
of Legendre polynomials, with the normalisation: 



the factor multiplying /(«)) can be expanded into a Taylor 
series at v = with expansion coefficients ac 

ryf v 3 f 3 1 (l-ji 1 vr+ 1 

R(V) = t / d vHvj^r— . , , 

w/ AJ J w 7 2 (1 - LLv) a + 2 

( 1_ 7 2 (1 -p.v)(l -fliv)) } 

L 

v 2 dv d/2i d^ c f(v) ^2 ai(a, v),Vi,v) 



M 

E 

i=0 



(A5) 



We choose the direction of the outgoing photon as z-axis. 
The cosine of the polar angle of the electron is then fa 
and the azimuth of the electron direction is denoted by 
4> c . With this choice of coordinate system, the integration 
over electron direction {dfa d<f> e ) becomes trivial. The re- 
maining electron velocity integral can be expressed in the 
form of moments, which can easily be calculated numeri- 
cally: 



p i (n)P j (n)d l i = 



Sn ■ 

2i + l 



(A2) 



These polynomials obey an addition theorem (see e.g. 
Landau & Lifschitz |1988| Eq. (c.10)): 



Pi(v) 



Pl(»)PM 



(A3) 



m)! 



{l + m)\ 



PTMPrfai) cos(m(^-^)). 



where (j), fi = cos 8 and 4>\ , fi± = cos 9\ define two direc- 
tions with 7 = arccos^ the angle between these directions. 
Equation (Al) can then be written: 



i f R(Ti)dt = ^w i (a,Q)Pi(rinfa 1 )=:K(j i ,fi 1 )< t k4) 

71 J i=0 

To calculate the coefficients uii(a,Q) for given values 
of the temperature 0, we have to perform the integral over 
the electron velocity v, which is tractable in the limit of 
high or low electron plasma temperature. In these limits 
we can e xpan d the integrand into a series in v, as described 
in Sect. A.l for t he n on-relativistic case, and in 1/7, as 
described in Sect. A. 2 for the relativistic case. 



The expansion of the angular and spatial dependence 
of J((i, t) due to Legendre and Chcbyshcv polynomials is 
described in Sects. A. 3 and A. 4. 



A.l. Expansion of the phase function in the limit O <g; 1 

For small plasma temperatures the main contribution to 
the integral over electron velocity arises from the region of 
small velocities. Therefore, the integrand (more precisely 



(v k ) := 4vr / f(v)v k+2 dv. 



(A6) 



The normalisation of f(v) gives (v°) — 1, and (v 2 ) = 39 
for non-relativistic electron plasma temperatures. We used 



Mathematica (see e.g. Wolfram 1991 ) to expand the inte- 
grand of R(rj) up to 16th. order in v. Therefore (v k ) has 
to be calculated for k = 0, 2, 4, . . . , 16. The odd moments 
vanish due to the angular integration. This expansion in 
v results in a polynomial series of order M = L/2 + 2 in 
77. (Note, however, that one may still choose to truncate 
at M < L/2 + 2). As an example we give the expansion 
up to 4th. order in v. The coefficients are then given by: 



wo(a,6) 



uj 2 (a, 9) 



u 3 (a, 9) 



1 



<^ 2 ) , 2 



150 
2 



(a 2 



+ 3a) 
3a) (7a 2 



21a + 22), 



(w 2 )(a 2 + 3a + l) 



-<^(2a4- 
25 V 

1 , (V 2 ) , 2 

2 + — {a 



12a 3 + 21a 2 + 9a + 6), 



+ 3a - 6) 



210 



(10a 4 + 60a 3 + 55a 2 - 105a + 78) . 



-L^(a-l)(a + 4) 



50 



(a~l)(a + 4)(a 2 + 3a-7) 



w 4 (a,9) = ^-f(a-l)(a-2)(a + 4)(a + 5). 
175 



(A7) 
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For very small electron temperatures and sufficiently small 
values of a, the sum of Eq. (A5) converges very quickly. 
If we truncate it at (e.g.) L = 2, the coefficient of the 
isotropic part of Eq. (A4) is 



w (a, 9) = l + 9(a 2 + 3a), 



(A8) 



in agreement with Eq. (18) of Titarchuk & Lyubarskij 
(|l995l)il. 



A. 2. Expansion of the phase function in the limit 0»1 

In the relativistic limit, the expansion has to be done in a 
slightly different way, because the phase function contains 
singular parts at v = 1. Using the orthogonality of the 
Legendre polynomials, the coefficients of Eq. (Al) can be 
written: 



Wi(a,6) 



2i+l 

47T 



(A9) 



where df2 — dr/dip with </> the azimuth of the incoming 
photon with respect to the outgoing one, which defines 
the z-axis. In this reference frame, the electron volume- 
element can be written d 3 v = v 2 dvdCl c . Using this rela- 
tion, and the definition of the phase function, the expan- 
sion coefficients can be written as 



w J (a,0) = 37r J dvv 



(A10) 



with the temperature independent kernel 



&»(a,7) 



2i + l 



Di\ a+2 J_ 
D ) Dt 

{i + (v') 2 }^(v) dndn c . 



(All) 



The axis of integration in this reference frame can be 
changed, so that the electron direction becomes the z-axis. 
This is expressed by df2df2 = dQdQi. With this choice 
of coordinates, it is easy to perform a Lorentz boost to the 
electron rest frame. The azimuths of the in- and outgoing 
photons do not change, and the angle between the photon 
direction and the boost direction changes according to 



t L = 



t' 



1 + vp> 

The differential solid angle transforms as 



dQ 



dil' 



7 2 (1 + vp') 2 



3 Equation (A. 16) 

<% = ^{l + yKa + 3) 



should 



(A12) 



(A13) 



read: 



6]} 



These transformations lead to 
2i + l 



Wi(a,7) 



(4tt) 2 7 2 



dp,' dp[ d<t> d^ 



where the argument of the Legendre polynomials has to 
be inserted as 



77 = 1- 



1 



1 



7 2 (l + vp')(l + vp{) 



(A15) 



To solve the above integrals one might use integral tables, 
or computer programs. We used Mathematica to find ana- 
lytic solutions for i = 0, 1, 2, 3 (see Titarchuk & Lyubarskij 
199i| E q s - (A10)-(A15), for i = 0). In the limit v -> 1 



these solutions diverge. Therefore we separated the lead- 
ing order terms. Extracting a 7-dependent factor accord- 
ing to: 

^(a, 7 ) = (2 7 ) 2Q+2 ^(a), 



these solutions are: 

a(a + 3) 



u>o(a) 
&i(a) 
0)2(0:) 



= -3 



= 5 a 



(a + l)(a + 2) 2 (a + 3) ' 
a(a + 3) + 4 



(A16) 
(A17) 



(a + 2) 2 (a + 3) 2 ' 

a(a + 3) + 4 
(a + 2) 2 (a + 3) 2 (a + 4) ' 

a(a + 3) 



w 3 (a) = 7a(l-a) ? — - 

(a + 2)^(a + 3)^(a 4- 4) (a + 5) 

Again, the isotropic part, (1)0(0,7), is in agreement with 



Eq. (A17) of Titarchuk & Lyubarskij (1995) 



For high plasma temperatures one can expand the inte- 
grand of Eq. (A10) into a series. Using the following rela- 
tion: 



^ d «=£v 1_ ? d7= £ (1 ~27 + -" )d7 (A18) 

Eq. (A10) becomes: 
3 1 



u>i(a) 



4 9K 2 (l/9) 

OO 



(A19) 



The above integration^ can be expressed in terms of the 
incomplete T-functionD. Then the coefficients Wj(a, 9) can 
be written: 

(29) 2q 

Ui(a,Q) = 3 o)i(a) 



K 2 (l/9) 

e 2 r(2c 



1 \ 1 / In 

3>e)-2 r ( 2a + 1 '9^ 



(A20) 
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Together with Eqs. (A17) this defines the expansion coef- 
ficients (for i = 0,1,2,3) of the phase function for a high vectors P. Let W be tjhe matrix, which transforms the 
electron plasma temperature. 



with a diagonal, M + 1-dimensional matrix to(a, 0) and 
vectors P. Let W be t 
Legendre polynomial as 

P(2/ii - 1) = W Pfa) 



A. 3. The angular dependence of J(/j,,t) 
The angular dependence of the intensity J(/x, r) can be Then ' for ^ - °> one can write 



expressed in the form of a series of Legendre polynomials 
with r-dependent coefficients. Because of the symmetry of 
the boundary conditions, it is necessary only to calculate 
r) for one half-space: < fi < 1 and —To < t < tq. 
The intensity in the second half-space is then given by 



(A21) 



Having in mind the orthogonality relation for Legendre 
polynomials on the interval (e.g.) < /i < 1 : 



P n (2/i-l)P fe (2^-l)dM = 



1 



2n+ 1 



0~nk , 



(A22) 



we expand the intensity as follows 

N 



E 

N 

E 



2n + 1 



P n (2/x + l)Q-(r) 



2n + 1 



P n (2/*-l)Q+(r). 



(A23) 



Defining the new variables £ and £, according to: 
C := 2/i + l| 



£ := 2/i - 1 



p<0 ' 
p>0 ' 



(A24) 



the inversion of Eq. (A23) can be written as 
Qn(r) = 



J(^,r)P n (C)dC : 



1 

ftM - / 



J 



)P 1 (0dC- 



(A25) 



Noting that P n {^) = (— l) n Pn( — /•*)> the symmetry condi- 
tion (Eq. A21) gives 



;(r) = (-irQ+(-r). 



if(^Mi) = P(M)w(a,e)P( m ) 



u(a,e)P0<i) = w(a,e)W- 1 WP(/n) 
= : w+(a,e)P(2^i - 1) 



(A28) 



(A29) 



with an M + l~dimensional, non-diagonal, matrix 
w + (a, O). An analogous transformation in the range /ii < 
leads to 

m) = P(A0[w + (a,e)P(2/ii-l)H(/*i) 

+ aT (a, 6) P(2/*i + 1) H(-^i) ] , (A30) 
with the Heaviside function: 



H(/xi) 



/*! > 

Mi < 



(A31) 



The source function B(h,t) (Eq. 9) becomes with this 
(and the use of Eq. (A25), and suppressing the depen- 
dence of w ± on a and O): 

P(/x,r) = \[P{^)uj + Q + (t) + P(jjl) u -Q-(t)] . (A32) 

Inserting this into the Boltzmann equation (Eq. 8), sub- 
stituting £ = 2/i — 1 for /i > 0, multiplying by Pfe(£), and 
integrating £ from -1 to 1, we obtain 



— 5— flbCO^l— S— >t) 



(t) + (ay Q, (r) 



i 

d£Pfc (£)•/( 



(A33) 



The indices i,j, k run from to N. If N > M the matrices 
oj'lj and w« can be defined with all elements equal to for 
i,j > M. UN < M only the elements with i,j < N 
are taken into account, so that in both cases the matrices 
(A26) formally become N + 1-dimensional. Using the recurrence 
relation for Legendre polynomials, 



It is useful to rewrite the kernel K(fi, //i) using a transfor- 
mation of variables and to split it into two parts according £ p k (£) 
to Hi < and /ii > 0. This is done in the following way: 
written in the form of a scalar product, the right-hand 
side of Eq. (A4) reads: 



L Pk , (0 k 



1 



2fc + 1 



2k + 1 



(A34) 



5 The matrix W can be calculated from Eq. (A28) by ex- 
pressing the vectors P(fti) and P(2/ii — 1) in the common 



basis (1, h, fx 2 , . 
(A27) ofP(Mi). 



^ ) and inverting the matrix of coefficients 
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the left hand side of Eq. (A33) can be written in terms of arises compared to the above treatment of the angular de- 
a tridiagonal matrix, whose elements are defined as: pendencc in that we have to take into account the bound- 

ary condition, which is that no radiation enter the disk 
Zu := S k ,i +1 + S k i + 8 k ,i-i • (A35) f ™ ° ut side: 

J(M,To)|^ <0 =0 = J^-ro)!^. (A42) 



We introduce at this point a new variable according to 



We further define a matrix with the elements: 

i T 
S kl := \ I ^ Pk ( 0Pi (i±l) . (A36) = - • (A43) 



TO 

_1 Using this transformation and Eq. (A23), the boundary 

TT . in - • -ci ;i n ri /aoo\ condition reads for all expansion coefficients: 

Using these definitions, and again Eq. (A25), Eq. (A33) 1 

becomes Q+(-l) = = Q~(l) . (A44) 

N BO + ( ) N ^ e cnose an expansion of Q + in Chebyshev polynomials 

Y, Z *-^p- = E 5 « ['"tjQti'r) + WyQjM] according to 

1=0 ' i,j=0 K 

-2Qt(r). (A37) Q+(y) = qo + J2[l-T i (-y)]q i . (A45) 

»=1 

Having in mind the symmetry relation, Eq. (A26), and ... ,~ , > , .... „ , knn \ 

, „ . ,. . .. '-iv n with q iy\ g lvcn by the symmetry condition Eq. (A26). 

denning a diagonal matrix according to , ', .... . _ T , . : 

00 ° lhe boundary condition then gives q = 0. Inserting this 

E jm := (-iy 6 jm , (A38) into Ec l' ( A41 ) one § ets 

Eq. (A37) reads in matrix form l^-i^+i ^M zft = ^ T^y) [(-l) i+1 M - M] qi 

z dQpr)_ = s r + +{t) _ +{ _ t) , k 

dr V yJ +V[M + M]q 4 . (A46) 

-2Q+(t). (A39) i= i 

This can be expressed in somewhat shorter way, using the Multiplying this equation by Tj (y)/y/l-y 2 and integrat- 

definitions in § f from ~ 1 to 1 > lt becomes 



m : = -2ii + sc+, Ly^i+i iKy ]°y m 7 qi = 

To £-f J WT - -y 2 

M := Sco-E. (A40) - 1 

x 1 

Now the Boltzmann equation can be written as f ^jiv) [7 i _ |\/|1 . 
7 9Q+(r) 



= MQ+(t) + MQ+(-t), (A41) 



^ + V (&y^jM=[M + M\ Qi . (A47) 

and represents a system of N+ 1 coupled differential equa- r— J J ^/l — y 2 

tions. _1 

For a given temperature in the non-relativistic regime This is a vector equation with if x K matrices, defined 

(9 < 1), the moments (v k ) are to be calculated up to k = by the integrals over Chebyshev polynomials. We denote 

L. The matrix elements of M and M are then polynomials these matrices by 

in a of maximal order L. l „ , . a „ . , 

/T(v)—T(v) 
3Ky 'dy l ^fJ_ 
\J\ — y 2 

function (see Eq. (A20)). 1 



l 



A. 4. The spatial dependence of J{h,t) ^ '~ J u 

Equation (A41) is a system of coupled differential equa- 



2 

1 
1 



tions for the N + 1 expansion coefficients of J (fx, t), which . . f Tj(y) 

contain the r-dependence of the intensity. In expanding ^ h '- 7 ' t ' J J f\ _ ~p. ' 
these coefficients themselves into a series, the problem is _1 

reduced to an algebraic eigenvalue problem. A difference (E t )ifc := (— Sik ■ (A48) 



12 



U.D.J. Gieseler & J.G. Kirk: An eigenfunction method for the comptonisation problem 



The matrix elements are straightforwardly calculated us- 
ing the orthogonality relation for Chebyshev polynomials: 



dy TMTM 

V 1 -y 2 




i = j = 



(A49) 



the expression for the derivative: 

dTjjy) = -iyTjjjD+iTj-^y) 
Ay (1 — y 2 ) 

and the recurrence relation 

T i+1 (y) = 2yT i (y)-T i _ 1 (y). 



(A50) 



(A51) 



The boundary condition qo = is not explicitly in- 
cluded in the above matrix equation, but must be inserted 
'by hand' into Eq. (A47). This is done by extending the 
matrix equation to i,j = 0, . . . , K, and adding the equa- 
tions 



K 



J2 C H U 

*=0 

with 

(Qn = 



o 



j = K and 
otherwise . 



(A52) 



(A53) 



All other matrices have to be set to for j — K. The 
matrices are then given by 



(Td). 



i > j and i + j 
otherwise , 



odd 



(To) 



(Th 



i ^ j or j = K 
i = j ^ and j ^ K 
i=j=0 , 

3=0 

otherwise . 



(A54) 



It is useful to write Eq. (A47) in a different way. It is 
a vector equation in which the elements are themselves 
vectors, multiplied by matrices. This form of matrix mul- 
tiplication is just the outer product of matrices, denoted 
by <S>. The K + 1 vectors of dimension N + 1 can be writ- 
ten as a common vector in a (K + 1) • (N + l)-dimensional 
product-space: 



WJ 



(A55) 



Note that the vector q does not have an index but in- 
cludes all expansion coefficients of the angular and spatial 



dependence of J(/x,r). Equation (A47) now reads 

[-(E t T d )®Z]q = [C®1 + T h ® (M + M) 
Tq — 

+ (E t T Q ) <g> M - T ®M]q.(A56) 

Let us denote the matrix on the left hand side, which does 
not depend on any parameter, by 



A := (E t T d )<8>Z, 



(A57) 



and the right hand matrix, which depends on a and 9 
(due tow ± (a,e) in M and M) by 0). Equation (A56) 
reads with these definitions: 

[— A - D(a,0)l q = F(a,0,T O )q = 0. (A58) 
To - — _ - 

I : defines an algebraic eigenvalue problem for all of the 
expansion coefficients q. 
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